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Partially Penalized Immersed Finite Element Methods 
for Parabolic Interface Problems 
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Abstract We present partially penalized immersed finite element methods for solving parabolic 
interface problems on Cartesian meshes. Typical semi-discrete and fully discrete schemes are 
discussed. Error estimates in an energy norm are derived. Numerical examples are provided to 
support theoretical analysis. 
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1 Introduction 


In this article, we consider the following parabolic equation with the Dirichlet boundary condition 

— ~V-(pVu) = f(X,t), X = (x,y) € fi+Ufi", t€ (0,T], (1.1) 

u\an = g(X,t), t G (0, T], (1.2) 

u\t=o = uq(X), X G O. (1.3) 


Here, 0 is a rectangular domain or a union of several rectangular domains in M 2 . The interface T C fi is 
a smooth curve separating Q into two sub-domains and such that H = Q~ U U T, see the left 
plot in Figure 1. The diffusion coefficient f3 is discontinuous across the interface, and it is assumed to be 
a piecewise constant function such that 




rr, xen-, 

\P + , xen+, 


and min{/3 , /3 + } > 0. We assume that the exact solution u to the above initial boundary value problem 
satisfies the following jump conditions across the interface T: 


Mr — 


Pp- 


= 0 . 


u r 


(1.4) 

(1.5) 
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Figure 1: The simulation domain fl (left), body-fitting mesh (middle), and non-body-fitting mesh (right). 


Interface problems appear in many applications of engineering and science; therefore, it is of great 
importance to solve interface problems efficiently. When conventional finite element methods are em¬ 
ployed to solve interface problems, body-fitting meshes (see the mid plot in Figure 1) have to be used 
in order to guarantee their optimal convergence [1, 2, 3, 5]. Such a restriction hinders their applications 
in some situations because it prevents the use of Cartesian mesh unless the interface has a very simple 
geometry such as an axis-parallel straight line. Recently, immersed finite element (IFE) methods have 
been developed to overcome such a limitation of traditional finite element methods for solving interface 
problems, see [6, 7, 8, 9, 12, 13, 14, 15, 18, 20, 23]. The main feature of IFE methods is that they 
can use interface independent meshes; hence, structured or even Cartesian meshes can be used to solve 
problems with nontrivial interface geometry (see the right plot in Figure 1). Most of IFE methods are 
developed for stationary interface problems. There are a few literatures of IFEs on time-dependent inter¬ 
face problems. For instance, an immersed Eulerian-Lagrangian localized adjoint method was developed 
to treat transient advection-diffusion equations with interfaces in [22], In [19], IFE methods were applied 
to parabolic interface problem together with the Laplacian transform. Parabolic problems with moving 
interfaces were considered in [10, 16, 17] where Crank-Nicolson-type fully discrete IFE methods and IFE 
method of lines were derived through the Galerkin formulation. 

For elliptic interface problems, classic IFE methods in Galerkin formulation [13, 14, 15] can usually 
converge to the exact solution with optimal order in H 1 and L 2 norm. Recently, the authors in [18, 23] 
observed that their orders of convergence in both H 1 and L 2 norms can sometimes deteriorate when the 
mesh size becomes very small, and this order degeneration might be the consequence of the discontinuity 
of IFE functions across interface edges (edges intersected with the interface). Note that IFE functions in 
[13, 14, 15] are constructed so that they are continuous within each interface element. On the boundary 
of an interface element, the continuity of these IFE functions is only imposed on two endpoints of each 
edge. This guarantees the continuity of IFE functions on non-interface edges. However, an IFE function 
is a piecewise polynomial on each interface edge; hence it is usually discontinuous on interface edges. 
This discontinuity depends on the interface location and the jump of coefficients, and could be large 
for certain configuration of interface element and diffusion coefficient. When the mesh is refined, the 
number of interface elements becomes larger, and such discontinuity over interface edges might be a 
factor negatively impacting on the global convergence. 

To overcome the order degeneration of convergence, a partially penalized immersed finite element 
(PPIFE) formulation was introduced in [18, 23]. In the new formulation, additional stabilization terms 
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generated on interface edges are added to the finite element equations that can penalize the discontinuity 
of IFE functions across interface edges. Since the number of interface edges is much smaller than the total 
number of elements of a Cartesian mesh, the computational cost for generating those partial penalty terms 
is negligible. For elliptic interface problems, the PPIFE methods can effectively reduce errors around 
interfaces; hence, maintain the optimal convergence rates under mesh refinement without degeneration. 

Our goal here is to develop PPIFE methods for the parabolic interface problem (1.1) - (1.5) and to 
derive the a priori error estimates for these methods. We present the semi-discrete method and two 
prototypical fully discrete methods, i.e., the backward Euler method and Crank-Nicolson method in 
Section 2. In Section 3, the a priori error estimates are derived for these methods which indicate the 
optimal convergence from the point of view of polynomials used in the involved IFE subspaces. Finally, 
numerical examples are provided in Section 4 to validate the theoretical estimates. 

In the discussion below, we will use a few general assumptions and notations. First, from now on, we 
will tacitly assume that the interface problem has a homogeneous boundary condition, i.e., g = 0 for the 
simplicity of presentation. The methods and related analysis can be easily extended to problems with 
a non-homogeneous boundary condition through a standard procedure. Second, we will adopt standard 
notations and norms of Sobolev spaces. For r > 1 , we define the following function spaces: 

H r {n) = {u : v\ n s € H r (n s ), s = + or -} 


equipped with the norm 


IMfer (n) = M ffr(n-) + IMIffr (n +)> w G 

For a function z(X,t ) with space variable X = (x,y) and time variable t, we consider it as a mapping 
from the time interval [0, T] to a normed space V equipped with the norm || • ||y. In particular, for an 
integer A; > 1, we define 


L k (0,T;V) = < z : [0, T] —>• V measurable, such that / ||z(- 


, t)\\ydt < oo 


with 

IMlL fc (0,T;V) = \\ z {'ltyWvdt 

Also, for V = H r (Q), we will use the standard function space H p (0, T; FT r (fi)) for p > 0, r > 1. 

In addition, we will use C with or without subscript to denote a generic positive constant which may 
have different values according to its occurrence. For simplicity, we will use ut, utt, etc., to denote the 
partial derivatives of a function u with respect to the time variable t. 



2 Partially Penalized Immersed Finite Element Methods 

In this section, we first derive a weak formulation of the parabolic interface problem (1.1) - (1.5) based 
on Cartesian meshes. Then we recall bilinear IFE functions and spaces defined on rectangular meshes 
from [8, 15]. The construction of linear IFE functions on triangular meshes is similar, so we refer to 
[13, 14] for more details. Finally, we introduce the partially penalized immersed finite element methods 
for the parabolic interface problem. 
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2.1 Weak Form on Continuous Level 

Let Tfi be a Cartesian (either triangular or rectangular) mesh consisting of elements whose diameters 
are not larger than h. We denote by Mh and £h the set of all vertices and edges in Th, respectively. 
The set of all interior edges are denoted by £/,. If an element is cut by the interface T, we call it an 
interface element; otherwise, it is called a non-interface element. Let 72 be the set of interface elements 
and 77 l be the set of non-interface elements. Similarly, we define the set of interface edges and the set of 
non-interface edges which are denoted by £ l h and £%, respectively. Also, we use £ l h and £JJ to denote the 
set of interior interface edges and interior non-interface edges, respectively. With the assumption that 
r C ft we have £\ = £ l h . 

We assign a unit normal vector ng to every edge B £ £\ l . If B is an interior edge, we let Kb .i and 
Kb, 2 be the two elements that share the common edge B and we assume that the normal vector n# is 
oriented from Kb,i to Kb, 2 - For a function u defined on Kb ,i U Kb, 2 , we set its average and jump on B 
as follows 

"Mb = 2 {( u \k b ,i)\b + (u\k B ' 2 )\b) , Mb = (u\k b ,i)\b ~ (u\k b , 2 )\b- 
If B is on the boundary 9ft, n b is taken to be the unit outward vector normal to 9ft, and we let 

{Mb = Mb = u \b- 

For simplicity, we often drop the subscript B from these notations if there is no danger to cause any 
confusions. 

Without loss of generality, we assume that the interface F intersects with the edge of each interface 
element K £ 72 at two points. We then partition K into two sub-elements K~ and K + by the line 
segment connecting these two interface points, see the illustration given in Figure 2. 

To describe a weak form for the parabolic interface problem, we introduce the following space 

Vh = {v £ L 2 (ft) : v satisfies conditions (HV1) - (HV4)} (2-1) 

(HV1) v\ K G H 1 (K),v\ K s £ H 2 {K s ),s = ±,\/K £ %. 

(HV2) v is continuous at every X £ Mh ■ 

(HV3) v is continuous across each B £ £%. 

(HV4) v\ dn = 0. 


Note that functions in V), are allowed to be discontinuous on interface edges. We now derive a weak 
form with the space V/, for the parabolic interface problem (1.1) - (1.5). First, we assume that its exact 
solution u is in i/ 2 (ft). Then, we multiply equation (1.1) by a test function v £ Vjj and integrate both 
sides on each element K £ Th- If K is a non-interface element, a direct application of Green’s formula 


leads to 


/ utvdX + / /3Xu ■ XvdX — 
Ik Jk 


/ /3Vu • n xvds = / fvdX. 

‘dK Jk 


If K is an interface element, we assume that the interface F intersects dK at points D and E. Then, 


without loss of generality, we assume that T and the line DE divide K into up to four sub-elements, see 
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Figure 2: An interface element 


the illustration in Figure 2 for a rectangle interface element, such that 

k = {n + n k + ) u (rr n k~) u (n + n k~) u (ft - n k + ) 

Now, applying Green’s formula separately on these four sub-elements, we get 

- [ V • (fiXu)vdX 
Jk 

= - [ V • (/ 3 + Xu)vdX - [ V • (, 0~Xu)vdX 

Jn+nK+ Jn-nK- 

- [ V • ( 0 + Xu)vdX - [ V • (, 0-Vu)vdX 

Jn+nK - Jn~nK+ 

f3 + Xu • XvdX - 


in+nK+ 


id{n+rK+) 


f3 + Xu ■ nvds + (3 Xu ■ XvdX — j3 Xu ■ nvds 

Jn-rK- Jd{n-nK~) 


+ / 0 + Xu-XvdX- 
Jn+nK- 


/ /3 + Xu-nvds + / f3 Xu ■ XvdX — 

'd(n+n K-) Jn~nK+ 


'd(n~n K+) 


[3 Xu ■ nvds 


= / (3Xu ■ XvdX — / f3Xu ■ nvds — 


IK 


l OK 


Ik nr 


IfiXu ■ n] A - nr vds 


= / (3Xu ■ XvdX — / j3Xu ■ nvds. 


IK 


/ dK 


(2.3) 


The last equality is due to the interface jump condition (1.5). The derivation of (2.3) implies that (2.2) 
also holds on interface elements. 


Remark 2.1. Figure 2 is a typical configuration of an interface element. If the interface is smooth 
enough and the mesh size is sufficiently small, an interface is usually divided into three sub-elements, 
i.e., one of the two terms n K~ and n K + is an empty set. In this case, the related discussion is 
similar but slightly simpler. 


Summarizing (2.2) over all elements indicates 

f utvdX + V"' f (3Xu ■ XvdX — f I^fJXu • ns]} [u] ds = f fvdX. 
Jn Jk Jb Jn 


(2.4) 


Best 
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Let Hh = H 2 (Q) + Vh on which we introduce a bilinear form ci t : H\ x Hh 


a e (w,v)= [ PVv-VivdX— Y'' f j[/3Vu; • n^]} [v] ds 

,■ Jk Jb 

K Besi 


+e XJ J ^ 3Vv ' M ds + XJ 

Ber h B Beil' 


a 


B 


M N ds, 


IB 


(2.5) 


where a > 0, <To > 0 and \B\ means the length of B. Note that the regularity of u leads to 


e ^2 I • n B | l'u] ds = 0, ^ / tw- M f«1 ds = 0. 


Be£i 


IB 


Best 


IB 


We also define the following linear form 


L{y) = f fvdX. 

Jn 

Finally, we have the following weak form of the parabolic interface problem (1.1)-(1.5): find u : [0, T] —> 
H 2 (Q) that satisfies (1.4), (1.5), and 

(u t ,v)+ a e (u,v) = L(v), Vu G 14, (2.6) 

u(X,0) = uo(X), VXen. (2.7) 


2.2 Immersed Finite Element Functions 

In this subsection, to be self-contained, we recall IFE spaces that approximate 14. We describe the 
bilinear IFE space with a little more details, and refer readers to [13, 14] for corresponding descriptions 
of the linear IFE space on a triangular Cartesian mesh. Since an IFE space uses standard finite element 
functions on each non-interface element, we will focus on the presentation of IFE functions on interface 
elements. 

The bilinear (Qi) immersed finite element functions were introduced in [8, 15]. On each interface 
element, a local IFE space uses IFE functions in the form of piecewise bilinear polynomials constructed 
according to interface jump conditions. Specifically, we partition each interface element K = OA 1 A 2 A 3 A 4 
into two sub-elements K~ and K + by the line connecting points D and E where the interface T intersects 
with dK, see Figure 3 for illustrations. Then we construct four bilinear IFE shape functions = 
1, 2, 3,4 associated with the vertices of K such that 

0, y) = <4 + bf x + cfy + djxy, 
y) = a T +b~x + c~ y + d~ xy, 

according to the following constraints: 

• nodal value condition: 

4>i(Aj) = Sij, i,j = 1,2,3,4. (2.9) 



if (x,y) G K+, 
if (x,y) € K~, 
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Figure 3: Type I and Type II interface rectangles 


continuity on DE 


ir 02 A. n 


[^(^)i = o, ium = o, 

continuity of normal component of flux 

( IT r)rh; H 

ds = 0. 


dxdy 


= 0. 


IDE 


on 


( 2 . 10 ) 


( 2 . 11 ) 


It has been shown [7, 8] that conditions specified in (2.9) - (2.11) can uniquely determine these shape 
functions. Figure 4 provides a comparison of FE and IFE shape functions. 



Figure 4: Bilinear FE/IFE local basis functions. From left: FE, IFE (Type I), IFE (Type II) 


Similarly, see more details in [13, 14], on a triangular interface element K = AA 1 A 2 A 3 , we can 
construct three linear IFE shape functions 4>i,i = 1,2,3 that satisfy the first two equations in (2.10), 
(2.11), and 


— $iji i,j — 1,2,3. 

These IFE shape functions possess a few notable properties such as their consistence with the corre¬ 
sponding standard Lagrange type FE shape functions and their formation of partition of unity. We refer 
readers to [7, 8, 13, 23] for more details. 
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Then, on each element K G Th, we define the local IFE space as follows: 


S h (K) = span{cj)i , 1 < i < d K }, d K 


3, if K is a triangular element, 

4, if K is a rectangular element, 


where fa, 1 < i < dK are the standard linear or bilinear Lagrange type FE shape functions for K G Tfa, 
otherwise, they are the IFE shape functions described above. Finally, the IFE spaces on the whole 
solution domain Q are defined as follows: 


S h (Q) = {v£V h : v\ K E S h (K), MK G T h }. 



Figure 5: Bilinear FE (left) and IFE (right) global basis functions 


Remark 2.2. We note that an IFE function may not be continuous across the element boundary that 
intersects with the interface. An IFE shape function is usually not zero on an interface edge, see the values 
on the edge between the points (0,0) and (0,1) for the two IFE shape functions plotted in Figure 4■ On this 
interface edge, the shape functions vanish at two endpoints, but not on the entire edge. The maximum of 
the absolute values of the shape on that edge is determined by the geometrical and material configuration 
on an interface element. When the local IFE shape functions are put together to form a Lagrange type 
global IFE basis function associated with a node in a mesh, it is inevitably to be discontinuous on interface 
edges in elements around that node, as illustrated by the cracks in a global IFE basis function plotted in 
Figure 5. As observed in [18, 23], this discontinuity on interface edges might be a factor causing the 
deterioration of the convergence of classic IFE solution around the interface, and this motivates us to 
add partial penalty on interface edges for alleviating this adversary. 

2.3 Partially Penalized Immersed Finite Element Methods 

In this subsection we use the global IFE space Sh(Ll) to discretize the weak form (2.6) and (2.7) for 
the parabolic interface problem. While the standard semi-discrete or many fully discrete frameworks can 
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be applied, we will focus on the following prototypical schemes because of their popularity. 

A semi-discrete PPIFE method: Find Uh '■ [0, T\ —> 5/(0) such that 

(uh,t, Vh ) + a e {u h ,v h ) = L(v h ), \/v h E S h (Q,), (2.12) 

n ft (A,0)=%(A), VAeH, (2.13) 

where Uh is an approximation of uq in the space Sffifil). According to the analysis to be carried out in 
the next section, Uh can be chosen as the interpolation of uq or the elliptic projection of uq in the IFE 
space >S)j(Q). 


A fully discrete PPIFE method: For a positive integer Nt, we let At = T/7V) which is the time step 
and let t n = nAt, for integer n > 0. Also, for a sequence ip n , n > 1, we let 


dw n = 


ip 


n— 1 


At 


Then, the fully discrete PPIFE method is to find a sequence { uf} ] of functions in Sh{D) such that 

(dtK,v h ) + a e {6u n h + (1 - 6 )ul~ 1 ,v h ) = 0L n {v h ) + (1 - 6 )L n ~\v h ), Vv h E S h (n), (2.14) 

u° h {X) = u h (X), VlGfi. (2.15) 


Here, L n (yh) = f n f(X,t n )vh(X)dX,n > 0 and 6 is a parameter chosen from [0,1]. Popular choices for 
9 are 9 = 0, 9 = 1 and 0 = 1/2 representing the forward Euler method, the backward Euler method, and 
the Crank-Nicolson method, respectively. 

Remark 2.3. The bilinear form a e (-,-) in (2.6) is almost the same as that used in the interior penalty 
DG finite element methods for the standard elliptic boundary value problem [4, 11, 21] except that it 
contains integrals over interface edges instead of all the edges. This is why we call IFE methods based on 
this bilinear form partially penalized IFE (PPIFE) methods. As suggested by DG finite element methods, 
the parameter e in this bilinear form is usually chosen as —1, 0, or 1. Note that a e (-,-) is symmetric if 
e = — 1 and is nonsymmetric otherwise. 


3 Error Estimations for PPIFE Methods 

The goal of this section is to derive the a priori error estimates for the PPIFE methods developed in 
the previous section. As usual, without loss of generality for error estimation, we assume that g(X, t) = 0 
in the boundary condition (1.2) and assume F n 3D = 0. The error bounds will be given in an energy 
norm that is equivalent to the standard semi -H l norm. These error bounds show that these PPIFE 
methods converge optimally with respect to the polynomials employed. 
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3.1 Some Preliminary Estimates 

First, for every » 6 14, we define its energy norm as follows: 



For an element K G Tj,, let \K\ denote the area of K. It is well known that the following trace inequalities 
[21] hold: 

Lemma 3.1. There exists a constant C independent of h such that for every K G Th, 

\\v\\ L *(B) < C|S| 1 / 2 |K|- 1/2 (||u|| l2W + h\\Vv\\ L 2 (K) ), Vu g H\K), B c dK, (3.1) 

||Vu|| L2(B) < C'|5| 1 / 2 |lF|- 1 / 2 ( ||Vu|| i;2( x ) + h\\V 2 v\\ LHK) ), \/v G H 2 {K), B C dK. (3.2) 

Since the local IFE space Sh{K ) C H l {K) for all K e Th (e.g. [7, 8, 13]), the trace inequality (3.1) 
is valid for all v G Sh(K). However, for I\ e 7)(, a function v G Sh{K) does not belong to H 2 (K ) in 
general. So the second trace inequality (3.2) cannot be directly applied to IFE functions. Nevertheless, 
for linear and bilinear IFE functions, the corresponding trace inequalities have been established in [18]. 
The related results are summarized in the following lemma. 

Lemma 3.2. There exists a constant C independent of interface location and h but depending on the 
ratio of coefficients (3 + and (3~ such that for every linear or bilinear IFE function v on K G Tf, 

WPvdh^B) < Ch 1 / 2 \K\- 1 / 2 \\'/f 3 X 7 v\\ L 2 ( K) , Vu G S h (K), B C dK, d = x or y, (3.3) 

II PVv • nfl|| l2(b) < Ch l / 2 \K\- l / 2 \\JpVv\\ L 2 {K) , Vu G S h (I<),B C dK. (3.4) 

As in [18], using Young’s inequality, trace inequalities and the definition of || • \\h, we can prove the 
coercivity of the bilinear form a e (-, •) on the IFE space Sh(Q) with respect to the energy norm || • ||^. The 
result is stated in the lemma below. 

Lemma 3.3. There exists a constant k > 0 such that 

ae(v h ,v h ) > K\\v h \\l, \/v h eS h (n) (3.5) 

holds for e = 1 unconditionally and holds for e = 0 or e = — 1 when the penalty parameter a Q B in a e (-, •) 
is large enough and a > 1. 

For any t G [0, T], we define the elliptic projection of the exact solution u(-,t ) as the IFE function 
Uh(-,t) G S h (Q) by 

a e (u - u h , v h ) = 0, Mv h G S h {Vl). (3.6) 

Lemma 3.4. Assume the exact solution u is in H 2 { 0, T ; H A ft)) and a = 1. Then there exists a constant 
C such that for every t G [0, T] the following error estimates hold 


II U ^/l||/l — Ch\\u\\ H 3^y 

(3.7) 

\\{u-Uh)t\\h < Ch\\u t \\ S 3 {n y 

(3.8) 

\\{u - u h ) tt \\ h < Ch\\u t t\\ S 3( Q y 

(3.9) 
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Proof. First, the estimate (3.7) follows directly from the estimate derived for the PPIFE methods for 
elliptic problems in [18]. Because of the linearity of the bilinear form, we have that 


a e ((u - u h ) t ,v h ) 



Uh,v h ) = 0, \/v h eS h (n). 


This indicates that the time derivative of the elliptic projection is the elliptic projection of the time 
derivative. Thus, for any given t € [0,T], ut € the estimate (3.8) follows from the estimate 

derived for the PPIFE methods for elliptic problems in [18] again. Similarly, we can obtain (3.9). □ 


3.2 Error estimation for the semi-discrete method 

The a priori error estimates for semi-discrete PPIFE method (2.12)-(2.13) for parabolic interface 
problem is given in the following theorem. 

Theorem 3.1. Assume that the exact solution u to the parabolic interface problem (1.1)-(1.5) is in 
i7 1 (0, T; H 3 (Q)) for e = —1 and in H 2 (0,T; H 3 (Q)) when e = 0, 1, and uq G H 3 {Fl). Let Uh be the 
PPIFE solution defined by semi-discrete method (2.12)-(2.13) with a = 1 and Uh(-, 0) = «/,,(•) being the 
elliptic projection of uq. Then there exists a constant C such that 

\[u(-,t) — Uh(-,t)\\h < Ch^\\uo\\jj 3 ^ + 11ut||£,2(o i T;h 3 (r2)))’ Vt > 0, (3.10) 

for e = —1, and 

IK*,*) - u h (-,t)\\h 

— + IK(*> 0)llil3(Q) + \\ u t\\ L2(o,T;ii 3 (n)) + \\ Uu \\l 2 (o,t-,h 3 (Q))) > — 0, (3-11) 

for e = 0 or 1. 

Proof. Let iih be the elliptic projection of u defined by (3.6) and we use it to split the error u — Uh into 
two terms: u — Uh = r) — £ with rj = u — Uh and £ = — uh- For the first term, by (3.7), we have the 

following estimate: 

IK*,*)IU<C7i|K*>*)lltf3(n) ^ Ch(\\u 0 \\ 

< C7i(lKI|/j-3( n ) + IIKlL 2 (o,T;i73(n)))- ( 3 -12) 

Then, we proceed to bound \\i\\h- From (2.6), (2.12) and (3.6), we can see that £ satisfies the following 
equation: 

(£t, v h ) + a e (£, v h ) = ( rj t , v h ), \/v h G S h (Q). (3.13) 

Choosing Vh = it in (3.13), we have 

||£ t || 2 + a e (£,£ t ) = (3.14) 

If e = — 1, using the symmetry property of a e (-, •), Cauchy-Schwarz inequality and Young’s inequality 
iu (3.14), we get 

1161,2 + - |W 11611 - c “’"ii 2 + \ 1161,2 ■ (3 - 15 > 
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For any t £ (0, T], integrating both sides of (3.15) from 0 to t. using the fact £(-,0) = 0 and (3.8), we 
obtain ^ ^ 

\M 2 + \aM-,tU(-,t))<cJ o ||?ft|| 2 dt < Ch 2 j |KII^3 (n) ■ ( 3 - 16 ) 

Using coercivity of a e (-, •) in (3.16), we have 

11 11L 2 (0,£;L 2 (f2)) + ||£||/i < Ch\\llt llL 2 (0,T;fl 3 (Q))' (3.17) 

Finally, applying the triangle inequality, (3.12) and (3.17) to u — Uh = rj — £ leads to (3.10). 

When e = 1 or 0, a e (-, •) is not symmetric. However, we have 

( 0 &) = \ 0 + 1 (o e (£, 6) - °e (6,0) 

>2^ ( ^0-c ,| l6 | U || e | U 

(3.i8) 

Substituting (3.18) into (3.14) and then integrating it from 0 to t, we can get 

7)j Ut\\ 2 dr + -n\\0\l<C J (\\Vt \\ 2 + Ut\\l + U\\h) dT - (3.19) 

Now we need the bound of HCtlU- From (3.13), we can easily get 

(Ot,Vh) + a e (£ t ,Vh) = (vtt,Vh), Vv h £S h (Q),t> 0. (3.20) 

Choose Vh = 0 in (3.20) and use the coercivity of a e (-, •) to get 

2^H&I| 2 + K ll&llft — + ll^ll 2 )- 

Integrating the above inequality from 0 to t and using the Gronwall inequality, we obtain 

T HtWldr < C f || vttfdr + C||6(-,0)|| 2 . (3.21) 

Jo Jo 

Let t = 0 and then choose Vh = £t(-,0) in (3.13) to get 

||6(-,0)||<|M-,0)||. (3.22) 

Substituting (3.21) and (3.22) into (3.19) and then using the Gronwall inequality again, we obtain 

I* \M 2 dr + U\\l < C [\\\vtW 2 + htt\\ 2 )dr + C\\r, t (; 0)|| 2 . 

Jo Jo 

Applying (3.8) and (3.9) to the above yields 

ll£i II.L 2 (0,£;L 2 (f2)) + ll£IU — Ch(\\ut(', 0)||^ 3 (Q)) + 11 Hi II L 2 (0,T;i7 3 (Q)) + \\ u tt H,L 2 (0,T;i7 3 (fi))) ' (3.23) 

Finally, applying the triangle inequality, (3.12) and (3.23) to u — Uh = g — f yields (3.11). □ 

Remark 3.1. By slightly modifying the proof for Theorem 3.1 we can show that estimates (3.10) and 
(3.11) still hold when Uh is chosen to be the IFE interpolation of uq. 
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3.3 Error estimation for fully discrete methods 

In all the discussion from now on, we assume that u9 = uh is the elliptic projection of uq in the initial 
condition for the parabolic interface problem. Also, for a function 4>(t), we let (f) n = cj)(t n ),n > 0. 


3.3.1 Backward Euler method 

The backward Euler method corresponds to the method described by (2.14) with 9 = 1. From (2.6), 
(2.14) and (3.6), we get 

(dtC, v h ) + a e {£ n , v h ) = ( d t if , v h ) + (r n , v h ), \/v h <E S h (Q), (3.24) 

where r n = — (it™ — dtu n ). We choose the test function Vh = dt£ n in (3.24) and use the Cauchy-Schwarz 
inequality on the right hand side to obtain 

m n \\ 2 + a e (C,dtO < (l \dtv n \\ + ||r n ||)||^n < (ll V|| 2 + |H| 2 ) + 111<9^'|| 2 - (3-25) 

There are three cases depending on the parameter e. We start from the case in which e = — 1. By the 
symmetry and the coercivity of the bilinear form a e (-, •), we have 


a«(r,9 t c)=^»<(r.c-c- 1 ) 


C 1 ) 


Thus, we have 

)ii9.n 2 + T i (a«(c,n-««K n - 1 ,r- 1) ) < iiai”ii 2 + 

Multiply (3.26) by 2A t and then sum over n to get 


n 11 2 i n r ni|2 


n= 1 


A t J2 m n f + ae(Z k , a k ) < 2At ^ ( \\d tV n \\ 2 + 

n= 1 

By Holder’s inequality and (3.8), we have 

! | 2 - 


II2 


, , n n — n n ~ l \ 2 

- n " 2 - 1 1 1 - ? - | dX = 

At 


Jn 

1 


rt n 


-f 

!i 'Ai J t n -1 

h 2 r 


dX 


- a* i-, I™ dT 5 C A( 

Applying Taylor formula and Holder’s inequality, we have 

ft" ' 2 


2 ~ dr 

HZ(Q) aT - 


|r n || 2 = 


(3.26) 


(3.27) 


(3.28) 


[ \u™ - d t u n \ 2 dX = f ^ - f (t-t n 1 )u t tdt dX<^~ f \\u tt \\ 2 dT. (3.29) 

Jtt Jfl 2At J t n -1 O 
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Substituting (3.28) and (3.29) into (3.27) and then using the coercivity of a e (-, •), we obtain 


+ \\£ k \\h ^ < ^(^' 2 |l' Ut llL2( 0iT;i 5-3(Q)) + (^ i ) 2 |l u tillL2( 0iT;L 2(Q))j. (3.30) 

n =1 

Finally, applying the triangle inequality, (3.13) and (3.30) to u k — u k = ij k — £ k yields 

ll' u ' ~ u h\\h 5; C(MH U o|| #3(Q) + \\ u t\\L 2 (0,T;H 3 (fl))) + ^t\\ u tt ||L 2 (0,T;L 2 (Q))) (3.31) 

for any integer k > 0. 

Now we turn to the cases where e = 0 or e = 1 that make the bilinear form in the PPIFE methods 
nonsymmetric. We start from 

1 


ae(C,dtC) = Xt ae(C,C-C- 1 ) 


1 


> 


2 At 

1 

2A t 
1 

2A t 


2A t 


(a e (r,n - a e (r -1 ,r _1 )) + ^ t a ^ n ,c - e* -1 ) 
(a e (r,n - ^(r -1 ,^" 1 )) + l(a t (e,d t o - a^c.c- 1 ) 
(a e (r,n - aer-sr- 1 )) - c(\\d t c\\i + + lie 


ae{c-c-\c- 1 ) 


n l|2 

/i 


Substituting it into (3.25) leads to 

iiian 2 + ^(«<(nn - 0 < (£»- i ,5”- 1 )) 


< iianf + imf+c 

Multiply (3.32) by 2A t and sum over n to obtain 


"t 


+ iims + ni(). 


(3.32) 


E At n 8 <f”ii 2 +«k *6 s Emkwi 5 + i'*i , )+<7 E+ceA fiini (3.33) 


71=1 


71=1 


71=1 


71=1 


In order to bound we first derive from (3.24) that 

71=1 

i 


— (dtC - d t C 1 ,v h ) + a e (d t C, Vh) = (dttif, v h ) + ( d t r n , v h ), \/v h G S h (tt). (3.34) 


At 

Let Vh = dtC 1 in (3.34) to get 

1 


2A t 

Then we can easily obtain 


(Il9.ni 2 - l|9,5”- 1 || 2 ) + 48,a\l < (l|9,«"H + l|9 1 r”||)||a,n. 


9,{ t || 2 + E^I|9.niS<CEA((||9„. f ”|| 2 + ||9,r”|| 2 ) + C||8 1 { 1 || 2 . (3.35) 


71=2 


71=2 
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Let n = 1 and Vh = dt ? = £VAi in (3.24), then we have 


ll^ 1 !! 2 + — a e {i\e) < (ll^ll + lk 1 ||)||^ 1 ||. 


Thus 


At 


ll^ 1 || 2 + < C'(||a t r7 1 | |2 + Hr- 1 !! 2 ). 

Applying this to (3.35) yields 

k k 

E < C E A *(ll««u”ll 2 + Il9<r"ll 2 ) + Cdia.,, 1 II 2 + Hr 1 !! 2 ). 

n=l n=2 

Inserting (3.36) into (3.33), then applying the Gronwall inequality, we obtain 
k k 

ii«‘iij<cE A ‘(ii8.>?“ii 2 + ini 2 )+c’E A *(iis<ffl“ii 2 + ii«<'-"ii 2 )+c’(iivii 2 + ii>' 1 i 


n= 1 


n =2 


We now estimate the last four terms in (3.37). It is easily to see that 

' ij n — 2 ?y n_1 + ?7 n “ 


ll^ n f = 


(At) 2 


dX 


L (TW L m(t " ■- (W il - t)e!i ) “ 


< 


1 


3At Jtn-2 
This inequality and (3.9) lead to 

k 


la\\ 2 dt. 


J2 At WdttV n \\ 2 < Ch2 \\ u tt\\ 2 L 2 {0tT .H 3m y 


Also, we have 


n =2 


0 - 2« n " 1 + n n " 2 

<9 t r = * 


At 


(At) 2 

r tn i f t;n i r *" -1 

= / u tttdt-j—r 2 / Uut{t n l -t) 2 dt + / u m {t-t n 2 ) 2 dt,. 

Jtn-1 (At)- J t n-1 \Aty J t u -2 

Then, the application of Holder’s inequality leads to 

^ At||^r n || 2 < (At) 2 ^ f \\u m \\ 2 dt+l f \\u t u\\ 2 dt + ^ f ^ \\u ttt \\ 2 dt\ 

n =2 n=2 4t"-i Jt n -‘ 2 ) 

— < ~'{At) 2 \\ u ttt\\ 2 L2(0 : T;L 2 (n))- 


(3.36) 


. (3.37) 


(3.38) 


(3.39) 
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As for the last two terms on the right hand side of (3.37), we have 


'tfdt < h 2 — / \\u t \\% 3 , n ,dt 


and, by (3.29), 




P II 2 = j a k 1 - dtu'fdX < K (A J o lk,l| 2 *) 


(3.40) 


(3.41) 


Now, substituting (3.28), (3.29) and (3.38)-(3.41) into (3.37), we obtain 


ii^isc ik 


'L 2 ( 0 ,T-H 3 (Q)) 


+ \Wtt\ 


L 2 (0.T;//3(n)) ^ At 


4" A , / \\ u t\\ £r3rn^dt h 


+C ^||Kl| 2 (0,T;L 2 (Q)) + Kti 11| 2 (0,T;L 2 (n)) + J 0 \\ U tt\\ 2 dt'j (At) 2 . 

Again, applying the estimate for the triangle inequality and (3.12) to u k — uf = r] k — f k , we obtain 


\u k -4\\ h 


+ a, / 11^*11 h 


— C ^II M o||^-3(q) + 11 ^ 11L 2 (0,T: ,H 3 (fi)) + W Utt W L 2 (0,T-,H 3 (n)) + ll U *ll#3(! 

+C (\\utt\\L 2 ( 0 ,T-L 2 (n)) + 11 11L 2 ( 0 ,T ;L 2 (f 2 )) + \\utt\\ 2 dt^j ^ At. 


Now let us summarize the analysis above for the backward Euler PPIFE method in the following theorem. 
Theorem 3.2. Assume that the exact solution u to the parabolic interface problem (1.1)-(1.5) is in 

~ ~ ( 'j Nt 

i7 2 (0, T; H 3 (Q)) n H 3 (0, T; L 2 (Q)) and uq € H 3 (Q). Let the sequence uN be the solution to the 

l J 71=0 

backward Euler PPIFE method (2.14)-(2.15). Then, we have the following estimates: 

(1) Ife = -1, then there exists a positive constant C independent of h and At such that 

o<«<v f ~~ U ^ h ~ + H^llL 2 (o,T;A3(n))) + a *IMIl 2 (o,T;L 2 (Q))) • (3-42) 

(2) If e = 0 or 1, then there exists a positive constant C independent of h and At such that 


max K - KWh 

0 <n<Nt 


11^011 fr3(n\ A 


i7 3 (fi) + ll^hIIL 2 (0,T ;H 3 (Q)) + \\ Utt \\L 2 (0,T;H 3 (rt)) + J 


\M% 3 in) dt ) )h 


+C ^||KIl 2 (0,T;L 2 (T2)) + 11 u ttt 11 L 2 (0,T;L 2 (H)) + J \\ u tt\fdt^j ^ At. 


(3.43) 
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3.3.2 Crank-Nicolson method 


Now we conduct the error analysis for the Crank-Nicolson method which corresponds to 9 = 1/2 in 
(2.14). From (2.6), (2.14) and (3.6), we have 

(dtC, v h ) + \ a e(C + £ n- \ v h ) = ( dtif , v h ) + (r”, v h ) + (r^, v h ), Mv h G 5^(11), (3.44) 

where 

1/2 - i W + «r‘), rj = -(uT 1 ' 2 - d t u"). 

Taking Vh = <9t£ n = (£ n — £ n_1 )/A t in (3.44) and applying the Cauchy-Schwarz inequality, we get 

N n || 2 + 2- r- 1 ) < (ll VII + ||r?|| + ||r?||) ||$H 

< c(|| + || r i || 2 + || r 2 || 2 ) + \m n t- (3.45) 


If e = —1, due to the symmetry of a £ (-, •), we can rewrite (3.45) as 

ii^rii 2 + 2 ^(«e(r,n-°, ( r- 1 ,r- 1) )<c , (ii^i |2 

Multiplying (3.46) by 2A t and summing over n, we have 

k 

K\\e\\l < <Cj2^(\\dtV n \\ 2 + Kll 2 

n= 1 


II2 




(3.46) 


(3.47) 


We note that (3.28) is still a valid estimation for ||<9t7y n || 2 ; hence, we proceed to estimate ||r ? 1 l || 2 and |||| 2 
From the Taylor formula and Holder’s inequality, we obtain 


I r n 112 _ 
Kill — 


«r 1/2 -)«+<-■' 2 

r t n ~ 1 / 2 


dX 


[ -,[ [ u m {t-t n 1 )dt+ ( u tt t(t n - t)dt\ dX 

Jn 4 \ J t n- 1 Jtn- 1/2 / 


<C(At) 3 / |Ktt|| 2 df, 
Jt 11 -- 1 


(3.48) 


and 


r*o II 2 = 


£ (<" 1/2 -diU n )) 2 dA 
r i / r in_1/2 

L (W [L Um(t - + i-va ““ (t ” “ *> 


— t) dt dX 


<C(At) 3 \\u m \\' 2 dt. 

Jt n - 1 


(3.49) 


17 



Using (3.28), (3.48) and (3.49) in (3.47) yields 


U k \\l < ^(^ 2 |MlL2(0,T;i73(Q)) + ( At ) 4 |l u '«t|lL 2 (0,T;L 2 (n)))' 

Finally, we obtain an estimate for u k — u\ by applying the above estimate for f k , the triangle inequality 
and (3.7) to the splitting u k — ti, = y k — f k , and we summarize the result in the following theorem. 


Theorem 3.3. Assume that the exact solution u to the parabolic interface problem (1.1)-(1.5) is in 

~ ~ ( h Nt 

i7 1 (0, T; 77 3 (14)) n i7 3 (0, T; L 2 (i7)) and uq G Assume the sequence uU is the solution to 

l J 71=0 

the PPIFE Crank-Nicolson method (2.14)-(2.15) with e = —1. Then, there exists a positive constant C 
independent of h and At such that 

0<n<Nt ^ U _ U h\\h < C* (MII u o||#3(q) + \\ut\\ l 2 (0,T;H 3 (£1))^ + (^0 II u ttt\\L 2 (0,T;L 2 (Q.))J ■ (3.50) 

Remark 3.2. The choice of e = — 1 for the PPIFE Crank-Nicolson method is very natural because the 
method inherits the symmetry from the interface problem and its algebraic system is easier to solve. On 
the other hand, even though the non-symmetric PPIFE Crank-Nicolson methods based on the other two 
choices of e = 0 and e = 1 also seem to work well as demonstrated by the numerical results in the next 
section, the asymmetry in their bilinear forms hinders the estimation of several key terms in the error 
analysis so that the related convergence still remains elusive. 

Remark 3.3. We can replace the bilinear form a e (-, •) with the one used in the standard interior penalty 
DC finite element methods to obtain corresponding DGIFE methods for the parabolic interface problems. 
Furthermore, the error estimation for PPIFE methods can also be readily extended to the corresponding 
DGIFE methods. However, as usual, these DGIFE methods have much more unknowns than the PPIFE 
counterparts; hence they are less favorable unless features in DG formulation are desired. 


4 Numerical Examples 


In this section, we present some numerical results to demonstrate features of PPIFE methods for parabolic 
interface problems. 

Let the solution domain be fi = (0,1) x (0,1) and the time interval be [0,1]. The interface curve T 
is chosen to be an ellipse centered at the point (xq, yo) with semi-radius a and b, whose parametric form 
can be written as 

(x = x 0 + acos(0), / 4 -jx 

\y = yo + bsm(e). 

In our numerical experiments, we choose xq = yo = 0, a = 7r/4, b = 7r/6, and 9 G [0, 7r/2] . The interface 
T separates P into two sub-domains P _ = {(x,y) : r(x,y) < 1} and P + = {(x,y) : r(x,y) > 1} where 


r(x,y) 




+ 


(y - yo ) 2 

b 2 
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The exact solution for the parabolic interface problem is chosen to be 

if (x,y) g 

e 1 , if (x,y)eQ + , 

where p = 5 and the diffusion coefficients (3^ vary in different numerical experiments. 

We use a family of Cartesian meshes {Th, h > 0}, and each mesh is formed by partitioning 17 into 
N s x N s congruent squares of size h = 1/N S for a set of values of integer N s . For fully discretized methods, 
we divide the time interval [0,1] into Nt subintervals uniformly with t n = nAt, n = 0,1, • • • , JVj, and 
At = 1 /Nt. Also, we have observed that the condition numbers of the matrices associated with the 
bilinear forms in these IFE methods is proportional to /r -2 , similar to that of the standard finite element 
method; therefore, usual solvers can be applied to efficiently solve the sparse linear system in these IFE 
methods. 

First, we consider the case in which the diffusion coefficient (f3~, (3 + ) = (1,10) representing a moderate 
discontinuity across the interface. Both nonsymmetric (e = 1) and symmetric (e = —1) PPIFE methods 
are employed to solve the parabolic interface problem. For penalty parameters, we choose = 1 for 
the nonsymmetric method and <j \g = 100 for the symmetric method, while a = 1 for both methods. 
Both backward Euler and Crank-Nicolson methods are employed and the time step is chosen as At = 2 h. 
Errors of nonsymmetric and symmetric PPIFE backward Euler methods in L°°, L 2 and semi-# 1 norms 
are listed in Table 1 and Table 2, respectively. Errors of nonsymmetric and symmetric PPIFE Crank- 
Nicolon methods are listed in Table 3 and Table 4, respectively. All errors are computed at the final time 
level, i.e. t = 1. 


h 

II ■ IIl°° 

rate 

• IIl 2 

rate 

• Iff 1 

rate 

1/10 

2.7866#—2 


8.2619#—2 


2.1079#—0 


1/20 

7.9371#—3 

1.8118 

2.0935#—2 

1.9805 

1.0659#—0 

0.9838 

1/40 

2.6530#—3 

1.5810 

5.3984#—3 

1.9553 

5.3875#—1 

0.9844 

1/80 

8.9636#—4 

1.5655 

1.4473#—3 

1.8991 

2.7065#—1 

0.9932 

1/160 

3.3405#—4 

1.4240 

4.1586#—4 

1.7992 

1.3567#—1 

0.9963 

1/320 

1.3871#—4 

1.2680 

1.3204#—4 

1.6551 

6.7927#—2 

0.9980 

1/640 

6.2344#—5 

1.1538 

4.7909#—5 

1.4626 

3.3986#—2 

0.9990 

1/1280 

2.9411#—5 

1.0839 

1.9763#—5 

1.2775 

1.6998#—2 

0.9996 


u{t,x,y) = 




A+ 


rP _Ai 

0 + + 0 ~ 


Table 1: Errors of nonsymmetric PPIFE backward Euler solutions with (3 =1, j3 + = 10 at time t = 1 

In Table 1 and Table 2, we note that errors in semi-# 1 norms for both nonsymmetric and symmetric 
PPIFE backward Euler methods demonstrate an optimal convergence rate 0(h) + 0(At), which confirms 
our error estimates (3.42) and (3.43). Also note that the order of convergence in L 2 norm approaches 1 as 
we perform uniform mesh refinement. This is consistent with our expectation of the order of convergence 
0(h 2 ) + O(At) in L 2 norm although such an error bound has not been established yet. Errors gauged in 
L°° norm indicate a first order convergence for backward Euler method. 

In Table 3 and Table 4, the convergence rate in semi-# 1 norm confirms our error estimate (3.50) for 
Crank-Nicolson method. Moreover, errors in L 2 norm is of second order convergence which agrees with 


19 




h 

II ■ IIl°° 

rate 

• IIl 2 

rate 

• \m 

rate 

1/10 

6.6821E-2 


8.1952E-2 


2.1051E-0 


1/20 

1.5332E—2 

2.1237 

2.1070E-2 

1.9596 

1.0654E-0 

0.9826 

1/40 

5.1586E-3 

1.5715 

5.4326E-3 

1.9554 

5.3876E-1 

0.9836 

1/80 

1.5387E-3 

1.7453 

1.4582E-3 

1.8974 

2.7067E-1 

0.9931 

1/160 

4.9034E-4 

1.6498 

4.1727E-4 

1.8052 

1.3567E-1 

0.9964 

1/320 

1.7632E-4 

1.4755 

1.3212E-4 

1.6591 

6.7927E-2 

0.9980 

1/640 

7.1949E-5 

1.2932 

4.7927E-5 

1.4630 

3.3986E-2 

0.9990 

1/1280 

3.1775E-5 

1.1791 

1.9763E-5 

1.2780 

1.6998E-2 

0.9996 


Table 2: Errors of symmetric PPIFE backward Euler solutions with (3 =1, f3 + = 10 at time t = 1 


h 

• L°° 

rate 

II ' IIl 2 

rate 

1 • Ilp 

rate 

1/10 

5.1829E-2 


9.3610E-2 


2.1106E-0 


1/20 

1.0369E-2 

2.3215 

2.2475E-2 

2.0583 

1.0658E-0 

0.9857 

1/40 

2.8024E-3 

1.8875 

5.6292E-3 

1.9973 

5.3870E-1 

0.9843 

1/80 

7.1649E-4 

1.9676 

1.4091E-3 

1.9982 

2.7063E-1 

0.9931 

1/160 

1.7881E-4 

2.0026 

3.5445E-4 

1.9911 

1.3566E-1 

0.9963 

1/320 

4.5518E-5 

1.9739 

8.8742E-5 

1.9979 

6.7926E-2 

0.9980 

1/640 

1.1447E-5 

1.9914 

2.2156E-5 

2.0019 

3.3986E-2 

0.9990 

1/1280 

2.8833E-6 

1.9892 

5.5375E-7 

2.0004 

1.6998E-2 

0.9996 


Table 3: Errors of nonsymmetric PPIFE Crank-Nicolson solutions with f3 =1, f3 + = 10 at time t = 1 

our anticipated convergence rate 0(h 2 ) + 0(At 2 ). Errors in L°° norm also seem to maintain an optimal 
second order convergence. 

Next, we consider a larger discontinuity in the diffusion coefficient by choosing ((3~ ,/3 + ) = (1,10000). 
The nonsymmetric PPIFE method is used for spatial discretization in the experiment. We choose the 
penalty parameter cr ( / = 1 again for this large discontinuity case, since the coercivity bound is valid for 
any positive a^. Table 5 and Table 6 contain errors in backward Euler and Crank-Nicolson methods, 
respectively. Again, we observe that errors in semi -H 1 norm have an optimal convergence rate through 
mesh refinement for both methods. The convergence rate in L 2 norm is second order for Crank-Nicolson 
and first order for backward Euler. For symmetric PPIFE methods, we have observed similar behavior 
to the nonsymmetric methods provided that the penalty parameter <7^ is large enough. 
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h 

II ■ Hi 00 

rate 

• IIl 2 

rate 

• \m 

rate 

1/10 

1.0310E-1 


9.2384E-2 


2.1112E-0 


1/20 

1.4252E-2 

2.8447 

2.2543E-2 

2.0350 

1.0650E-0 

0.9872 

1/40 

4.2963E-3 

1.7401 

5.6546E-3 

1.9952 

5.3862E-1 

0.9836 

1/80 

1.0893E-4 

1.9796 

1.4190E-3 

1.9946 

2.7062E-1 

0.9930 

1/160 

2.8178E-4 

1.9508 

3.5605E-4 

1.9947 

1.3566E-1 

0.9963 

1/320 

6.6903E-5 

2.0744 

8.9021E-5 

1.9999 

6.7924E-2 

0.9980 

1/640 

1.6838E—5 

1.9903 

2.2251E-5 

2.0003 

3.3985E-2 

0.9990 

1/1280 

4.1878E-6 

2.0075 

5.5633E-6 

1.9999 

1.6998E-2 

0.9996 


Table 4: Errors of symmetric PPIFE Crank-Nicolson solutions with /3 =1, (3 + = 10 at time t = 1 


h 

• L°° 

rate 

' L 2 

rate 

1 ' I# 1 

rate 

1/10 

1.4637E-1 


4.7718E-2 


1.1268E-0 


1/20 

6.4974E-2 

1.1717 

1.6100E-2 

1.5675 

5.9288E-1 

0.9265 

1/40 

2.2137E-2 

1.5534 

4.3284E-3 

1.8951 

3.0548E-1 

0.9567 

1/80 

7.2728E-3 

1.6059 

8.4067E-4 

2.3642 

1.5187E-1 

1.0083 

1/160 

2.3746E-3 

1.6148 

2.0844E-4 

2.0119 

7.5576E-2 

1.0068 

1/320 

1.0006E-3 

1.2468 

5.2912E-5 

1.9779 

3.7807E-2 

0.9993 

1/640 

1.7030E-4 

2.5547 

1.4993E-5 

1.8193 

1.8900E-2 

1.0002 

1/1280 

6.3452E-5 

1.4244 

4.9410E-6 

1.6014 

9.4461E-3 

1.0006 
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